Description of two models:
Two simple models are implemented in ACT-R, performing the gambling task depending on two distinctive mechanisms.
The model randomly evaluates an action GUESS-MORE or GUESS-LESS, then retrieving memories about prior feedback. If “WIN” memory is retrieved, the model chooses the evaluated action. If “LOSE”/“Neutral” memory is retrieved, the model chooses the alternative action.
The model makes decision relying on procedural memory, pressing “MORE”/“LESS” key. GUESS-MORE and GUESS-LESS productions are competing. In the end, a feedback is given, the model receives either +1/-1/0 reward to all previous productions, and the utility of all previous fired productions will be affected.
model1.raw <- read.csv("./MODEL120210424_100307_fnca_log.csv") #%>% rename(Trial=X)
model2.raw <- read.csv("./MODEL220210424_100307_fnca_log.csv") #%>% rename(Trial=X)
#model1 <- read.csv("MODEL120210420.csv")
#model2 <- read.csv("MODEL220210420.csv")
model1.raw <- model1.raw %>%
mutate(CurrentResponse = case_when(Response=="j" ~ 3,
Response=="f" ~ 2),
RT = as.numeric(RT),
Epoch = as.numeric(Epoch))
model2.raw <- model2.raw %>%
mutate(CurrentResponse = case_when(Response=="j" ~ 3,
Response=="f" ~ 2),
RT = as.numeric(RT),
Epoch = as.numeric(Epoch))
# only
#model1 <- model1.raw %>% filter(ans==.1 & bll==.1 & lf==.1)
#model2 <- model2.raw %>% filter(alpha==.1 & egs==.1 & r==.1)
Win-stay probability(PSwitch) is calculated by: 1) count the number of ResponseSwitch grouped by BlockType and TrialType; 2) then calculate the proportion of ResponseSwitch (PSwitch) under groups.
future_moves <- function(responses) {
c(responses[2:length(responses)], NA)
}
past_moves <- function(responses) {
c(NA, responses[1:length(responses)-1])
}
count_responses <- function(model.dat) {
model.dat <- model.dat %>%
mutate(FutureResponse = future_moves(CurrentResponse),
PastResponse = past_moves(CurrentResponse),
PreviousFeedback = past_moves(TrialType),
ResponseSwitch = if_else(CurrentResponse == FutureResponse, 0, 1))
return(model.dat)
}
# NA for all first trial in each block
clean_previous_future <- function(model.dat) {
model.dat <- model.dat %>%
mutate(PastResponse=ifelse(BlockTrial==0, NA, PastResponse),
PreviousFeedback=ifelse(BlockTrial==0, NA, PreviousFeedback),
FutureResponse=ifelse(BlockTrial==7, NA, FutureResponse),
ResponseSwitch=ifelse(BlockTrial==7, NA, ResponseSwitch),
)
return(model.dat)
}
#model1.count <- clean_previous_future(count_responses(model1.raw))
Aggregate PSwitch, RT by ParamID, TrialType, BlockType
### This aggregate looks at overall, what is the average pattern of 50 simulation (50 epoch) across different parameter set
# aggregate data by parameter set, BlockType, TrialType
m1.aggregate.ByParam.CurrentResponse.BlockType <- function(m1.dat) {
res <- m1.dat %>%
filter(!is.na(ResponseSwitch) & !is.na(CurrentResponse)) %>%
group_by(ans, bll, lf, BlockType, TrialType) %>%
dplyr::summarise(PSwitch = mean(ResponseSwitch), RT = mean(RT)) %>%
ungroup() %>%
mutate(ParamID = group_indices(., ans, bll, lf))
return(res)
}
m2.aggregate.ByParam.CurrentResponse.BlockType <- function(m2.dat) {
res <- m2.dat %>%
filter(!is.na(ResponseSwitch) & !is.na(CurrentResponse)) %>%
group_by(egs, alpha, r, BlockType, TrialType) %>%
dplyr::summarise(PSwitch = mean(ResponseSwitch), RT = mean(RT)) %>%
ungroup() %>%
mutate(ParamID = group_indices(., egs, alpha, r))
return(res)
}
# aggregate data by parameter set, BlockType, PreviousFeedback
m1.aggregate.ByParam.PreviousFeedback.BlockType <- function(m1.dat) {
res <- m1.dat %>%
filter(!is.na(ResponseSwitch) & !is.na(PreviousFeedback) & !is.na(CurrentResponse)) %>%
group_by(ans, bll, lf, BlockType, PreviousFeedback) %>%
dplyr::summarise(PSwitch = mean(ResponseSwitch), RT = mean(RT))%>%
ungroup() %>%
mutate(ParamID = group_indices(., ans, bll, lf))
return(res)
}
m2.aggregate.ByParam.PreviousFeedback.BlockType <- function(m2.dat) {
res <- m2.dat %>%
filter(!is.na(ResponseSwitch) & !is.na(PreviousFeedback) & !is.na(CurrentResponse)) %>%
group_by(egs, alpha, r, BlockType, PreviousFeedback) %>%
dplyr::summarise(PSwitch = mean(ResponseSwitch), RT = mean(RT))%>%
ungroup() %>%
mutate(ParamID = group_indices(., egs, alpha, r))
return(res)
}
#model1.aggregate <- m1.aggregate.ByParam.CurrentResponse.BlockType(model1.count)
Aggregate PSwitch, RT by Epoch, TrialType, BlockType. (This only fit for one param set)
### This aggregate looks at for single parameter set, how consistent 50 runs are
# aggregate data by Epoch, BlockType, PreviousFeedback
aggregate.ByEpoch.CurrentResponse.BlockType <- function(m.dat) {
res <- m.dat %>%
filter(!is.na(ResponseSwitch) & !is.na(CurrentResponse)) %>%
group_by(Epoch, BlockType, TrialType) %>%
dplyr::summarise(PSwitch = mean(ResponseSwitch), RT = mean(RT))
return (res)
}
aggregate.ByEpoch.PreviousFeedback.BlockType <- function(m.dat) {
res <- m.dat %>%
filter(!is.na(ResponseSwitch) & !is.na(PreviousFeedback) & !is.na(CurrentResponse)) %>%
group_by(Epoch, BlockType, PreviousFeedback) %>%
dplyr::summarise(PSwitch = mean(ResponseSwitch), RT = mean(RT))
return (res)
}
#model1.aggregate <- m1.aggregate.ByParam.CurrentResponse.BlockType(model1.count)
Plot functions
# plot subject aggregated data, PSwitch by TrialType x BlockType
subj.PSwitch.aggplot <- function(subj.dat, subjID) {
res <- ggplot(subj.dat %>% filter(HCPID==subjID),
aes(x = TrialType, y = PSwitch)) +
facet_grid( ~ BlockType) +
geom_line(aes(group = BlockType), alpha=.5) +
geom_point(alpha=.5, size=3, aes(col=TrialType)) +
ylim(0,1) +
theme_pander()
return (res)
}
# plot model aggregated data, PSwitch by TrialType x BlockType x ParamID/Epoch
m.PSwitch.aggplot <- function(m.dat, byParam) {
if (byParam=="ParamID") {
res <- ggplot(m.dat, aes(x = TrialType, y = PSwitch, col = ParamID)) +
facet_grid( ~ BlockType) +
geom_line(aes(group = ParamID), alpha=.5) +
geom_point(alpha=.5, size=3) +
ylim(0,1) +
stat_summary(fun.data = "mean_cl_boot", col="black") +
theme_pander()
}
if (byParam=="Epoch") {
res <- ggplot(m.dat, aes(x = TrialType, y = PSwitch, col = Epoch)) +
facet_grid( ~ BlockType) +
geom_line(aes(group = Epoch), alpha=.5) +
geom_point(alpha=.5, size=3) +
ylim(0,1) +
stat_summary(fun.data = "mean_cl_boot", col="black") +
theme_pander()
}
return (res)
}
# plot model aggregated data, T by TrialType x BlockType x ParamID/Epoch
m.RT.aggplot <- function(m.dat, byParam) {
if (byParam=="ParamID") {
res <- ggplot(m.dat, aes(x = TrialType, y = RT, col = ParamID)) +
facet_grid( ~ BlockType) +
geom_line(aes(group = ParamID), alpha=.5) +
geom_point(alpha=.5, size=3) +
ylim(0,3500) +
stat_summary(fun.data = "mean_cl_boot", col="black") +
theme_pander()
}
if (byParam=="Epoch") {
res <- ggplot(m.dat, aes(x = TrialType, y = RT, col = Epoch)) +
facet_grid( ~ BlockType) +
geom_line(aes(group = Epoch), alpha=.5) +
geom_point(alpha=.5, size=3) +
ylim(0,3500) +
stat_summary(fun.data = "mean_cl_boot", col="black") +
theme_pander()
}
return (res)
}
#m.PSwitch.aggplot(model1.aggregate, "ParamID")
#m.RT.aggplot(model1.aggregate, "ParamID")
model1.count <- clean_previous_future(count_responses(model1.raw))
model2.count <- clean_previous_future(count_responses(model2.raw))
model1.aggregate <- m1.aggregate.ByParam.CurrentResponse.BlockType(model1.count)
model2.aggregate <- m2.aggregate.ByParam.CurrentResponse.BlockType(model2.count)
The following plots show the mean PSwitch and mean RT as a function of TrialType, BlockType. Each blue line represents a parameter set.
m.PSwitch.aggplot(model1.aggregate, "ParamID") + ggtitle("MODEL1")
## [1] "Choose one: ParamID or Epoch"
m.PSwitch.aggplot(model2.aggregate, "ParamID") + ggtitle("MODEL2")
## [1] "Choose one: ParamID or Epoch"
m.RT.aggplot(model1.aggregate, "ParamID") + ggtitle("MODEL1")
## [1] "Choose one: ParamID or Epoch"
m.RT.aggplot(model2.aggregate, "ParamID") + ggtitle("MODEL1")
## [1] "Choose one: ParamID or Epoch"
Model1:
Model2:
Model1:
Model2:
model1.pf_aggregate <- m1.aggregate.ByParam.PreviousFeedback.BlockType(model1.count)
model2.pf_aggregate <- m2.aggregate.ByParam.PreviousFeedback.BlockType(model2.count)
# plot subject aggregated data, PSwitch by PreviousFeedback x BlockType
subj.PSwitch.pfaggplot <- function(subj.dat, subjID) {
res <- ggplot(subj.dat %>% filter(HCPID==subjID),
aes(x = PreviousFeedback, y = PSwitch)) +
facet_grid( ~ BlockType) +
geom_line(aes(group = BlockType), alpha=.5) +
geom_point(alpha=.5, size=3, aes(col=PreviousFeedback)) +
ylim(0,1) +
theme_pander()
return (res)
}
# plot model aggregated data, PSwitch by PreviousFeedback x BlockType x ParamID/Epoch
m.PSwitch.pfaggplot <- function(m.dat, byParam) {
if (byParam=="ParamID") {
res <- ggplot(m.dat, aes(x = PreviousFeedback, y = PSwitch, col = ParamID)) +
facet_grid( ~ BlockType) +
geom_line(aes(group = ParamID), alpha=.5) +
geom_point(alpha=.5, size=3) +
ylim(0,1) +
stat_summary(fun.data = "mean_cl_boot", col="black") +
theme_pander()
}
if (byParam=="Epoch") {
res <- ggplot(model1.aggregate, aes(x = PreviousFeedback, y = PSwitch, col = Epoch)) +
facet_grid( ~ BlockType) +
geom_line(aes(group = Epoch), alpha=.5) +
geom_point(alpha=.5, size=3) +
ylim(0,1) +
stat_summary(fun.data = "mean_cl_boot", col="black") +
theme_pander()
} else {
print("Choose one: ParamID or Epoch")
}
return (res)
}
# plot model aggregated data, T by TrialType x BlockType x ParamID/Epoch
m.RT.pfaggplot <- function(m.dat, byParam) {
if (byParam=="ParamID") {
res <- ggplot(m.dat, aes(x = PreviousFeedback, y = RT, col = ParamID)) +
facet_grid( ~ BlockType) +
geom_line(aes(group = ParamID), alpha=.5) +
geom_point(alpha=.5, size=3) +
ylim(0,3500) +
stat_summary(fun.data = "mean_cl_boot", col="black") +
theme_pander()
}
if (byParam=="Epoch") {
res <- ggplot(model1.aggregate, aes(x = PreviousFeedback, y = RT, col = Epoch)) +
facet_grid( ~ BlockType) +
geom_line(aes(group = Epoch), alpha=.5) +
geom_point(alpha=.5, size=3) +
ylim(0,3500) +
stat_summary(fun.data = "mean_cl_boot", col="black") +
theme_pander()
} else {
print("Choose one: ParamID or Epoch")
}
return (res)
}
The following plots show the mean PSwitch and mean RT as a function of PreviousFeedback, BlockType. Each blue line represents a parameter set.
m.PSwitch.pfaggplot(model1.pf_aggregate, "ParamID") + ggtitle("MODEL1")
## [1] "Choose one: ParamID or Epoch"
m.PSwitch.pfaggplot(model2.pf_aggregate, "ParamID") + ggtitle("MODEL2")
## [1] "Choose one: ParamID or Epoch"
aov(PSwitch ~ BlockType * PreviousFeedback, model2.pf_aggregate %>% filter(PreviousFeedback!="Neutral")) %>%
anova_summary() %>%
xtable() %>%
kable(digits=4) %>%
kable_styling(bootstrap_options = c("striped", "hover"))
aov(RT ~ BlockType * PreviousFeedback, model2.pf_aggregate %>% filter(PreviousFeedback!="Neutral")) %>%
anova_summary() %>%
xtable() %>%
kable(digits=4) %>%
kable_styling(bootstrap_options = c("striped", "hover"))
Model1:
Model2:
Model1:
Model2:
blockId <- sort(rep(1:4, 8))
gambling <- gambling %>%
group_by(HCPID, RunNumber) %>%
mutate(BlockId = blockId,
Response = if_else(is.na(QuestionMark.RESP), 0, QuestionMark.RESP))
bg <- gambling %>%
group_by(BlockId, HCPID, RunNumber) %>%
summarise(BlockType = paste("Mostly", Mode(TrialType)))
gambling <- inner_join(gambling, bg)
past_moves <- function(responses) {
c(NA, responses[1:length(responses)-1])
}
gambling <- gambling %>%
group_by(RunNumber, BlockId) %>%
mutate(CurrentResponse = Response,
FutureResponse = future_moves(Response),
PastResponse = past_moves(Response),
PreviousFeedback = past_moves(TrialType),
RT = QuestionMark.RT/1200) %>%
filter(FutureResponse != 0) %>%
#filter(CurrentResponse != 0) %>%
mutate(ResponseSwitch = if_else(CurrentResponse == FutureResponse, 0, 1)) %>%
select(HCPID, RunNumber, BlockId, RT, Response, PastResponse,
CurrentResponse, FutureResponse, ResponseSwitch, BlockType, TrialType,
PreviousFeedback, FeedbackImage)
gambling$PreviousFeedback <- as_factor(gambling$PreviousFeedback)
gambling$PreviousFeedback <- revalue(gambling$PreviousFeedback,
c("1"="Reward",
"2"="Punishment",
"3"="Neutral"))
subj1.aggregate <- aggregate.CurrentTrial(gambling)
subj1.pf_aggregate <- aggregate.PreviousTrial(gambling)
PSwitch.aggplot(subj1.aggregate) + ggtitle("SUBJ - 100307")
RT.aggplot(subj1.aggregate) + ggtitle("SUBJ - 100307")
PSwitch.pfaggplot(subj1.pf_aggregate) + ggtitle("SUBJ - 100307")
RT.pfaggplot(subj1.pf_aggregate) + ggtitle("SUBJ - 100307")
gambling <- read_tsv("../../gambling_data.txt")
gambling.trials <- gambling %>%
#filter(RunNumber==1) %>%
select(HCPID, Trial, RunNumber, TrialType, RunTrialNumber, Block,
QuestionMark.RESP, QuestionMark.ACC, QuestionMark.RT,
R1MostlyReward1, R1MostlyReward2, R1MostlyPunishment1, R1MostlyPunishment3, QuestionMark.OnsetDelay,
ConsecSameResp, ConsecLargerGuesses, ConsecSmallerGuesses, ConsecRTLess200) %>%
rename(CurrentResponse=QuestionMark.RESP, RT=QuestionMark.RT) %>%
mutate(FutureResponse = future_moves(CurrentResponse),
PastResponse = past_moves(CurrentResponse),
ResponseSwitch = if_else(CurrentResponse == FutureResponse, 0, 1),
PreviousFeedback = past_moves(as.character(TrialType)),
BlockType = case_when(
!is.na(R1MostlyReward1) ~ "MostlyReward",
!is.na(R1MostlyReward2) ~ "MostlyReward",
!is.na(R1MostlyPunishment1) ~ "MostlyPunishment",
!is.na(R1MostlyPunishment3) ~ "MostlyPunishment"
)) %>%
select(-R1MostlyReward1, -R1MostlyReward2, -R1MostlyPunishment1, -R1MostlyPunishment3) %>%
filter(!is.na(TrialType))
#write.csv(gambling.trials, "../../gambling_clean_data.csv")
gambling.aggregate <- gambling.trials.session1 %>%
filter(!is.na(ResponseSwitch) & !is.na(RT)) %>%
group_by(HCPID, TrialType, BlockType, PreviousFeedback) %>%
summarise(PSwitch = mean(ResponseSwitch),RT = mean(RT))
ggplot(gambling.aggregate,
aes(x = TrialType, y = PSwitch, col = TrialType)) +
facet_grid( ~ BlockType, labeller = label_both) +
geom_point(position = position_jitter(width = 0.1, height = 0.05),
alpha = 0.1) +
scale_color_brewer(palette = "Set2") +
stat_summary(fun.data = "mean_cl_boot", col="black") +
theme_pander()
gambling.trials <- gambling.trials.session1 %>% select(HCPID, Trial, RunTrialNumber, TrialType, BlockType)
gambling.trials.session1 %>%
gghistogram("RT", add = "median")
gambling.clean <- read_csv("../../bin/gambling_clean_data.csv") %>%
select(-BlockType) %>% rename(BlockType=BlockTypeCoded)
gambling.cfagg <- gambling.clean %>%
na.omit() %>%
group_by(HCPID, BlockType, TrialType) %>%
summarise(PSwitch = mean(ResponseSwitch), RT = mean(RT))
gambling.pfagg <- gambling.clean %>%
na.omit() %>%
group_by(HCPID, BlockType, PreviousFeedback) %>%
summarise(PSwitch = mean(ResponseSwitch), RT = mean(RT))
# long format
gambling.cfagg <- left_join(gambling.cfagg %>% select(-RT) %>%
spread(key = TrialType, value = PSwitch, sep = "."),
gambling.cfagg %>% select(-PSwitch) %>%
spread(key = TrialType, value = RT, sep = "."),
by = c("HCPID", "BlockType"), suffix=c(".PSwitch",".RT"))
gambling.pfagg <- left_join(gambling.pfagg %>% select(-RT) %>%
spread(key = PreviousFeedback, value = PSwitch, sep = "."),
gambling.pfagg %>% select(-PSwitch) %>%
spread(key = PreviousFeedback, value = RT, sep = "."),
by = c("HCPID", "BlockType"), suffix=c(".PSwitch",".RT"))
#write_csv(gambling.cfagg, "../../bin/gambling_cfagg.csv")
#write_csv(gambling.pfagg, "../../bin/gambling_pfagg.csv")
gambling.clean <- read_csv("../../bin/gambling_clean_data.csv") %>% mutate(BlockType=BlockTypeCoded)
subjID = "100307_fnca"
# Pswitch as a function of current trial type and block type
subj.cfaggregate = gambling.clean %>%
filter(!is.na(ResponseSwitch) & !is.na(CurrentResponse)) %>%
group_by(HCPID, BlockType, TrialType) %>%
dplyr::summarise(PSwitch = mean(ResponseSwitch, rm.na=T), RT = mean(RT, rm.na=T))
# Pswitch as a function of previous trial type and block type
subj.pfaggregate = gambling.clean %>%
filter(!is.na(ResponseSwitch) & !is.na(PreviousFeedback)) %>%
group_by(HCPID, BlockType, PreviousFeedback) %>%
dplyr::summarise(PSwitch = mean(ResponseSwitch, rm.na=T), RT = mean(RT, rm.na=T))
Load cleaned subject data and select one subj data = 100307_fnca Note: here, ResponseSwitch is 1 if: CurrentResponse != FutureResponse; is 0 if CurrentResponse==FutureResponse
# test whether prop is calculate correctly
subj.cfaggregate %>% filter(HCPID==subjID)
gambling.clean %>% filter(HCPID==subjID) %>%
filter(!is.na(ResponseSwitch)) %>%
group_by(BlockType, TrialType , ResponseSwitch) %>%
dplyr::summarise(n=n()) %>%
mutate(freq = n / sum(n)) %>%
xtable() %>% kable(digits=4) %>% kable_styling(bootstrap_options = c("striped", "hover"))
#>> yes correct
The following plots show subject Pswitch pattern as a function of current TrialType and BlockType.
subj.PSwitch.aggplot(subj.cfaggregate, subjID) + ggtitle(paste("SUB", subjID, sep = "-"))
subj.PSwitch.pfaggplot(subj.pfaggregate, subjID) + ggtitle(paste("SUB", subjID, sep = "-"))
Distribution of model1 parameters :ans, :bll, :lf
Note: This is not grid-search parameter fitting. The parameter space was chosen by optimization function by adding null_penalty (count of null responses)
model1.full.count <- clean_previous_future(count_responses(model1.raw))
model1.full.agg <- model1.full.count %>%
filter(!is.na(ResponseSwitch) & !is.na(CurrentResponse)) %>%
group_by(ans, bll, lf, BlockType, TrialType) %>%
dplyr::summarise(PSwitch = mean(ResponseSwitch), RT = mean(RT))
model1.full.agg %>% gghistogram(x="ans", binwidth = .05) + ggtitle("MODEL1: :ans distrirbutioin")
model1.full.agg %>% gghistogram(x="bll", binwidth = .05) + ggtitle("MODEL1: :bll distrirbutioin")
model1.full.agg %>% gghistogram(x="lf", binwidth = .05) + ggtitle("MODEL1: :lf distrirbutioin")
Distribution model2 parameter :egs, :alpha, and :r
dim(model2.raw)
## [1] 339200 11
model2.full.count <- clean_previous_future(count_responses(model2.raw))
model2.full.agg <- model2.full.count %>%
filter(!is.na(ResponseSwitch) & !is.na(CurrentResponse)) %>%
group_by(egs, alpha, r, BlockType, TrialType) %>%
dplyr::summarise(PSwitch = mean(ResponseSwitch), RT = mean(RT))
model2.full.agg %>% gghistogram(x="egs", binwidth = .05) + ggtitle("MODEL1: :egs distrirbutioin")
model2.full.agg %>% gghistogram(x="alpha", binwidth = .05) + ggtitle("MODEL1: :alpha distrirbutioin")
model2.full.agg %>% gghistogram(x="r", binwidth = .05) + ggtitle("MODEL1: :r distrirbutioin")
Next, we calalculate RMSE for each parameter combination, find the best parameter with min RMSE
find_best_parameter <- function(model, model.agg, sub.agg) {
if (model=="model1") {
m1.param <- unique(model.agg[,c('ans','bll','lf')])
#m1.param=m1.param[1:10,]
m1.param$RMSE = NA
m1.param$RMSE = as.numeric(m1.param$RMSE)
for (i in 1:nrow(m1.param)) {
m1.agg <- model.agg %>%
filter(ans==m1.param[[i,'ans']] & bll==m1.param[[i,'bll']] & lf==m1.param[[i,'lf']])
rmse_value <- rmse(m1.agg$PSwitch, sub.agg$PSwitch)
m1.param[[i, 'RMSE']] <- rmse_value
}
return (m1.param)
} else {
m2.param <- unique(model.agg[,c('egs','alpha','r')])
#m2.param=m2.param[1:10,]
m2.param$RMSE = NA
m2.param$RMSE = as.numeric(m2.param$RMSE)
for (i in 1:nrow(m2.param)) {
m2.agg <- model.agg %>%
filter(egs==m2.param[[i,'egs']] & alpha==m2.param[[i,'alpha']] & r==m2.param[[i,'r']])
rmse_value <- rmse(m2.agg$PSwitch, sub.agg$PSwitch)
m2.param[[i, 'RMSE']] <- rmse_value
}
return (m2.param)
}
}
m1.bestp <- find_best_parameter("model1", model1.full.agg, subj.cfaggregate)
m2.bestp <- find_best_parameter("model2", model2.full.agg, subj.cfaggregate)
m1.bestp$model = "model1"
m2.bestp$model = "model2"
m1.bestp <- m1.bestp[order(m1.bestp$RMSE),]
m2.bestp <- m2.bestp[order(m2.bestp$RMSE),]
m1.bestp %>% head() %>% xtable() %>% kable(digits=4) %>% kable_styling(bootstrap_options = c("striped", "hover"))
| ans | bll | lf | RMSE | model |
|---|---|---|---|---|
| 1.2265 | 0.3263 | 0.757 | 0.2738 | model1 |
| 1.2257 | 0.1000 | 0.100 | 0.2744 | model1 |
| 1.2380 | 0.1000 | 0.100 | 0.2753 | model1 |
| 1.2265 | 0.3263 | 0.739 | 0.2754 | model1 |
| 3.0902 | 0.1000 | 0.100 | 0.2758 | model1 |
| 1.2265 | 0.3263 | 0.729 | 0.2770 | model1 |
m2.bestp %>% head() %>% xtable() %>% kable(digits=4) %>% kable_styling(bootstrap_options = c("striped", "hover"))
| egs | alpha | r | RMSE | model |
|---|---|---|---|---|
| 0.7295 | 0.1038 | 0.6972 | 0.2730 | model2 |
| 0.7295 | 0.1038 | 2.5030 | 0.2730 | model2 |
| 2.1291 | 0.1122 | 5.4124 | 0.2732 | model2 |
| 2.1287 | 0.1122 | 7.5503 | 0.2738 | model2 |
| 0.7295 | 0.1038 | 1.4486 | 0.2739 | model2 |
| 0.7295 | 0.1044 | 0.1000 | 0.2740 | model2 |
By looking at min RMSE:
Best parameter chosen for model1 :ans=1.2264712 :bll=0.3262827, :lf=0.757
Best parameter chosen for model2 is :egs=0.7294902 ; alpha=0.1037721, :r=0.6972462
Having the optimal parameter set, we can run simulation 800 times again by setting optimal parameter to the model. The plot shows the pattern of 800 simulation runs (Note: only single parameter set was loaded here)
subj1 = "./MODEL220210425_100307_fnca_best800.csv"
subj2 = "./seed/MODEL220210425_100307_fnca_best800.csv"
model2.best800 <- read.csv(subj1) %>%
mutate(CurrentResponse = case_when(Response=="j" ~ 3,
Response=="f" ~ 2),
RT = as.numeric(RT),
Epoch = as.numeric(Epoch))
model2.best800.count <- clean_previous_future(count_responses(model2.best800))
model2.best800.aggregate <- aggregate.ByEpoch.CurrentResponse.BlockType(model2.best800.count)
m.PSwitch.aggplot(model2.best800.aggregate, "Epoch") + ggtitle("MODEL2: optimal parameter 800 runs",
subtitle = paste("SUB", subjID, sep="-"))
ggarrange(m.PSwitch.aggplot(model2.best800.aggregate %>% filter(Epoch==0), "Epoch"),
m.PSwitch.aggplot(model2.best800.aggregate %>% filter(Epoch==1), "Epoch"),
m.PSwitch.aggplot(model2.best800.aggregate %>% filter(Epoch==2), "Epoch"),
m.PSwitch.aggplot(model2.best800.aggregate %>% filter(Epoch==3), "Epoch"),
m.PSwitch.aggplot(model2.best800.aggregate %>% filter(Epoch==4), "Epoch"),
m.PSwitch.aggplot(model2.best800.aggregate %>% filter(Epoch==5), "Epoch"))
Next, the following plots show the effect of random parameter (:ANS and :EGS) on model behaviors.
m1.gsfiles = list.files(path = "./gs_check",
pattern = paste(paste("^MODEL1.*", subjID, sep=""), "*_gs.csv$", sep=""), full.names = T)
m2.gsfiles = list.files(path = "./gs_check",
pattern = paste(paste("^MODEL2.*", subjID, sep=""), "*_gs.csv$", sep=""), full.names = T)
m1.gs <- read.csv(m1.gsfiles[1]) %>%
mutate(CurrentResponse = case_when(Response=="j" ~ 3,
Response=="f" ~ 2), RT = as.numeric(RT),
Epoch = as.numeric(Epoch))
m2.gs <- read.csv(m2.gsfiles[1]) %>%
mutate(CurrentResponse = case_when(Response=="j" ~ 3,
Response=="f" ~ 2), RT = as.numeric(RT),
Epoch = as.numeric(Epoch))
m1.gs <- m1.aggregate.ByParam.CurrentResponse.BlockType(clean_previous_future(count_responses(m1.gs))) %>%
mutate(ans = as.factor(ans), ParamID = as.factor(ParamID))
m2.gs <- m2.aggregate.ByParam.CurrentResponse.BlockType(clean_previous_future(count_responses(m2.gs))) %>%
mutate(egs = as.factor(egs), ParamID = as.factor(ParamID))
ggplot(m1.gs, aes(x = TrialType, y = PSwitch, col = ans)) +
facet_grid( ~ BlockType) +
geom_line(aes(group = ans), alpha=.5) +
geom_point(alpha=.5, size=3) +
ylim(0,1) +
stat_summary(fun.data = "mean_cl_boot", col="black") +
theme_pander() +
ggtitle("MODEL1: :ans effect", subtitle = paste("sub", subjID, sep="-"))
ggplot(m2.gs, aes(x = TrialType, y = PSwitch, col = egs)) +
facet_grid( ~ BlockType) +
geom_line(aes(group = egs), alpha=.5) +
geom_point(alpha=.5, size=3) +
ylim(0,1) +
stat_summary(fun.data = "mean_cl_boot", col="black") +
theme_pander() +
ggtitle("MODEL2: :egs effect", subtitle = paste("sub", subjID, sep="-"))
The following analysis is to calculate the Log-Likelihood of P(m|x) given human data; the RMSE of model data and human data; The BIC scores calculated from both RSS and LL.
model2.best800.aggregate <- model2.best800.aggregate %>%
group_by(BlockType, TrialType) %>%
dplyr::summarise(PSwitch.m = mean(PSwitch), PSwitch.sd = sd(PSwitch))
subj1.aggregate <- subj.cfaggregate %>% filter(HCPID == subjID) %>% rename(PSwitch.subj = PSwitch)
# bic
bic.ll <- function(k, n, logL) {
return(k * log(n) - 2 * logL)
}
#= n log(RSS/n) + (p + 1) log(n)
bic.rmse <- function(k, n, rss) {
return(n * log(rss/n) - (k + 1) * log(n))
}
# calculate log-liklihood
model2.best800.LL <- left_join(model2.best800.aggregate, subj1.aggregate) %>%
mutate(PSwitch.z = (PSwitch.m-PSwitch.subj)/(PSwitch.sd),
PSwitch.probz = dnorm(PSwitch.z),
PSwitch.logprobz = log(PSwitch.probz)) %>%
group_by(HCPID) %>%
dplyr::summarise(PSwitch.LL = sum(PSwitch.logprobz), RMSE = rmse(PSwitch.m, PSwitch.subj), RSS = RMSE^2 * 6)
model2.best800.LL %>% mutate(BIC.LL = bic.ll(k=3, n=6, PSwitch.LL), BIC.rmse = bic.rmse(k=3, n=6, RSS)) %>%
xtable() %>%
kable(digits=4) %>%
kable_styling(bootstrap_options = c("striped", "hover"))
| HCPID | PSwitch.LL | RMSE | RSS | BIC.LL | BIC.rmse |
|---|---|---|---|---|---|
| 100307_fnca | -8.1648 | 0.1824 | 0.1996 | 21.7049 | -27.5868 |